Мотивация: метод моментов

Метод моментов — это ещё один способ, наряду с методом максимального правдоподобия, оценки параметров распределения по данным $x_1,\ldots,x_N$. Суть его в том, что мы выражаем через параметры распределения теоретические значения моментов $\mu_k = \mathbb{E}x^k$ нашей случайной величины, затем считаем их выборочные оценки $\widehat{\mu}_k = \frac1N\sum_ix_i^k$, приравниваем их все друг к другу и, решая полученную систему, находим оценки параметров.

Можно доказать, что полученные оценки являются состоятельными, хотя могут быть смещены.

Пример 1. Оценим параметры нормального распределения $\mathcal{N}(\mu, \sigma^2)$ с помощью метода моментов.

Теоретические моменты равны

$$\mu_1 = \mu,\quad\mu_2 = \sigma^2 + \mu^2$$

Запишем систему:

$$\begin{cases} \widehat{\mu} = \frac1N\sum_i x_i,\\ \widehat{\sigma}^2 + \widehat{\mu}^2 = \frac1N\sum_ix_i^2 \end{cases}$$

Из неё очевидным образом находим

$$\widehat{\mu} = \frac1N\sum_ix_i$$

$$\widehat{\sigma}^2 = \frac1N\sum_ix_i^2 - \left(\frac1N\sum_i x_i\right)^2=$$

$$=\frac1N\sum_i\left(x_i - \widehat{\mu}\right)^2$$

Легко видеть, что полученные оценки совпадают с оценками максимального правдоподобия

Пример 2. Оценим параметр $\mu$ логнормального распределения

$$p(x) = \frac1{x\sqrt{2\pi\sigma^2}}\exp\left(-\frac{(\log{x} - \mu)^2}{2\sigma^2}\right)$$

при известном $\sigma^2$. Будет ли оценка совпадать с оценкой, полученной с помощью метода максимального правдоподобия?

Теоретическое математическое ожидание равно $\exp\left(\mu + \frac{\sigma^2}2\right)$, откуда мы сразу находим оценку $\widehat{\mu} = \log\left(\frac1N\sum_i x_i\right) - \frac{\sigma^2}2$.

Теперь запишем логарифм правдоподобия:

$$l(X) = -\sum_i\log{x_i} - \sum_i\frac{(\log{x_i} - \mu)^2}{2\sigma^2} + const$$

Дифференцируя по $\mu$ и приравнивая производную к нулю, получаем

$$\widehat{\mu}_{MLE} = \frac1N\sum_i\log{x_i}$$

что вовсе не совпадает с оценкой выше.

Несколько приукрасив ситуацию, можно сделать вывод, что первые два выборочных момента позволяют если не править миром, то уверенно восстанавливать параметры распределений. А теперь давайте представим, что мы посчитали $\frac1N\sum_ix_i$ и $\frac1N\sum_ix_i^2$, а семейство распределений пока не выбрали.

Как же совершить этот судьбоносный выбор? Давайте посмотрим на следующие три семейства и подумаем, в каком из них мы бы стали искать распределение, зная его истинные матожидание и дисперсию?

Почему-то хочется сказать, что в первом. Почему? Второе не симметрично — но почему мы так думаем? Если мы выберем третье, то добавим дополнительную информацию как минимум о том, что у распределения конечный носитель. А с чего бы? У нас такой инфомации вроде бы нет.

Общая идея такова: мы будем искать распределение, которое удовлетворяет только явно заданным нами ограничениям и не отражает никакого дополнительного знания о нём. Но чтобы эти нестрогие рассуждения превратить в формулы, придётся немного обогатить наш математический аппарат и научиться измерять количество информации.

Энтропия и дивергенция Кульбака-Лейблера

Измерять «знание» можно с помощью энтропии Шэннона. Она определяется как

$$\color{#348FEA}{H(P) = -\sum_xP(x)\log{P(x)}}$$

для дискретного распределения и

$$\color{#348FEA}{H(p) = -\int p(x)\log{p(x)}dx}$$

для непрерывного. В классическом определении логарифм двоичный, хотя, конечно, варианты с разным основанием отличаются лишь умножением на константу.

Неформально можно представлять, что энтропия показывает, насколько сложно предсказать значение случайной величины. Чуть более строго — сколько в среднем бит нужно потратить, чтобы передать информацию о её значении.

Пример 1. Рассмотрим схему Бернулли с вероятностью успеха $p$. Энтропия её результата равна

$$-(1 - p)\cdot\log_2(1 - p) - p\cdot\log_2{p}$$

Давайте посмотрим на график этой функции:

Минимальное значение (нулевое) энтропия принимает при $p\in\{0,1\}$. В самом деле, для такого эксперимента мы всегда можем наверняка сказать, каков будет его исход; обращаясь к другой интерпретации — чтобы сообщить кому-то о результате эксперимента, достаточно $0$ бит (ведь получатель сообщения и так понимает, что вышло).

Максимальное значение принимается в точке $\frac12$, что вполне соответствует тому, что при $p=\frac12$ предсказать исход эксперимента сложнее всего.

Дополнение для ценителей математики.

Попробуем для этого простого примера объяснить, почему среднее число бит, необходимых для передачи информации об исходе эксперимента, выражается формулой с логарифмами.

Теперь пусть $p$ произвольно. Рассмотрим $N»1$ независимых испытаний $x_1,\ldots, x_N$; среди них будет $n_0\approx (1-p)N$ неудачных и $n_1\approx pN$ удачных. Посчитаем, сколько бит потребуется, чтобы закодировать последовательность $x_i$ для известных $n_0$ и $n_1$. Общее число таких последовательностей равно $C_N^{n_1} = \frac{N!}{n_0!n_1!}$, а чтобы закодировать каждую достаточно будет $\log_2\left(\frac{N!}{n_0!n_1!}\right)$ бит — это количество информации, содержащееся во всей последовательности. Таким образом, в среднем, чтобы закодировать результат одного испытания необходимо

$$\frac1N\log_2\left(\frac{N!}{n_0!n_1!}\right)$$

бит информации. Перепишем это выражение, использовав формулу Стирлинга $\log{N!}\approx N\log{N} - N$:

$$\frac1N\left(\log_2{N!} - \log_2{n_0!} - \log_2{n_1!}\right) \approx$$

$$\approx const\cdot\frac1N\left(N\log_2{N} - N - n_0\log_2{n_0} + n_0 - n_1\log_2{n_1} + n_1\right) =$$

$$=const\cdot\frac1N\left(N\log_2{N} - (1-p)N\log_2{(1-p)N} - pN\log_2{pN}\right) =$$

$$=const\cdot\left(-(1-p)\cdot\log_2(1-p) - p\log_2{p}\right)$$

Вот мы и вывели формулу энтропии!

Пример 2. Энтропия нормального распределения $\mathcal{N}(\mu, \sigma^2)$ равна $\frac12\log(2\pi\sigma^2) + \frac12$, и чем меньше дисперсия, тем меньше энтропия, что и логично: ведь когда дисперсия мала, значения сосредоточены возле матожидания, и они становятся менее «разнообразными».

Энтропия тесно связана с другим важным понятием из теории информации — дивергенцией Кульбака-Лейблера. Она определяется для $p(x)$ $q(x)$ как

$$\color{#348FEA}{KL(p\vert\vert q) = \int p(x)\log{\frac{p(x)}{q(x)}}dx}$$

в непрерывном случае и точно так же, но только с суммой вместо интеграла в дискретном.

Дивергенцию можно представить в виде разности:

$$KL(p\vert\vert q) = (-\int p(x)\log{q(x)}dx) - (-\int p(x)\log{p(x)}dx)$$

Вычитаемое — это энтропия, которая, как мы уже поняли, показывает, сколько в среднем бит требуется, чтобы закодировать значение случайной величины. Уменьшаемое похоже по виду, и можно показать, что оно говорит о том, сколько в среднем бит потребуется на кодирование случайной величины с плотностью $p$ алгоритмом, оптимизированным для кодирования случайной величины $q$.

Иными словами, дивергенция Кульбака-Лейблера говорит о том, насколько увеличится средняя длина кодов для значений $p$, если при настройке алгоритма кодирования вместо $p$ использовать $q$. Более подробно вы можете почитать, например, в этом посте.

Дивергенция Кульбака-Лейблера в некотором роде играет роль расстояния между распределениями. В частности, $KL(p\vert\vert q)\geqslant0$, причём дивергенция равна нулю, только если распределения совпадают почти всюду. Но при этом она не является симметричной: вообще говоря, $KL(p\vert\vert q)\ne KL(q\vert\vert p)$.

Вопрос на подумать. Пусть $p(x)$ — распределение, заданное на отрезке $[a, b]$. Выразите энтропию через дивергенцию Кульбака-Лейблера $p(x)$ с равномерным на отрезке распределением $q_U(x)=\frac1{b-a}\mathbb{I}_{[a,b]}(x)$.

Ответ (не открывайте сразу; сначала подумайте сами!)

Распишем дивергенцию:

$$KL(p\vert\vert q_U) = -\left(-\int_a^b p(x)\log{p(x)}dx\right) - \int_a^b p(x)\log{\underbrace{q(x)}_{=\frac1{b-a}\text{ на }[a,b]}}dx=$$

$$=\log(b-a) - H(p)$$

Аналогичное соотношение можно выписать и для распределения, заданного на конечном множестве.

Принцип максимальной энтропии

Теперь наконец мы готовы сформулировать, какие распределения мы хотим искать.

Принцип максимальной энтропии. Среди всех распределений на заданном носителе $\mathbb{X}$, удовлетворяющих условиям $\mathbb{E}u_1(x) = \mu_1$, ..., $\mathbb{E}u_k(x) = \mu_k$, где $u_i$ — некоторые функции, мы хотим иметь дело с тем, которое имеет наибольшую энтропию.

В самом деле, энтропия выражает нашу меру незнания о том, как ведёт себя распределение, и чем она больше — тем более «произвольное распределение», по крайней мере, в теории.

Давайте рассмотрим несколько примеров, которые помогут ещё лучше понять, почему некоторые распределения так популярны:

Пример 1. На конечном множестве $1,\ldots,n$ наибольшую энтропию имеет равномерное распределение (носитель — конечное множество из $n$ элементов, других ограничений нет).

Доказательство: Пусть $p_i$, $i=1,\ldots,n$ — некоторое распределение, $q_i=\frac1n$ — равномерное. Запишем их дивергенцию Кульбака-Лейблера:

$$KL(p\vert\vert q) = \sum_i p_i\log{p_i} - \sum_i p_i\log{q_i} =$$

$$= -H(p) + \log{n}\underbrace{\sum_ip_i}_{=1}$$

Так как дивергенция Кульбака-Лейблера всегда неотрицательна, получаем, что $H(p)\leqslant\log{n}$. При этом равенство возможно, только если распределения совпадают.

Пример 2. Среди распределений, заданных на всей вещественной прямой и имеющих заданные матожидание $\mu$ и дисперсию $\sigma^2$ наибольшую энтропию имеет нормальное распределение $\mathcal{N}(\mu,\sigma^2)$.

Доказательство: Пусть $p(x)$ — некоторое распределение, $q(x)\sim\mathcal{N}(\mu, \sigma^2)$. Запишем их дивергенцию Кульбака-Лейблера:

$$KL(p\vert\vert q) = \int p(x)\log{p(x)}dx - \int p(x)\log{q(x)}dx =$$

$$= -H(p) - \int p(x)\left(-\frac12\log(2\pi\sigma^2) - \frac1{2\sigma^2}(x - \mu)^2\right)dx =$$

$$= - H(p) +\frac12\log(2\pi\sigma^2)\cdot\underbrace{\int p(x)dx}_{=1} + \frac1{2\sigma^2}\underbrace{\int(x - \mu)^2p(x)dx}_{=\mathbb{V}p=\sigma^2} =$$

$$= - H(p) + \underbrace{\frac12\log(2\pi\sigma^2) + \frac12}_{=H(q)}$$

Так как дивергенция Кульбака-Лейблера всегда неотрицательна, получаем, что $H(p)\leqslant H(q)$. При этом равенство возможно, только если распределения $p$ и $q$ совпадают почти всюду, а с точки зрения теории вероятностей такие распределения различать не имеет смысла.

Пример 3. Среди распределений, заданных на множестве положительных вещественных чисел и имеющих заданное матожидание $\lambda$ наибольшую энтропию имеет показательное распределение с параметром $\frac1{\lambda}$ (его плотность равна $p(x) = \frac1{\lambda}\exp\left(-\frac1{\lambda}x\right)\mathbb{I}_{(0;+\infty)}(x)$).

Все хорошо знакомые нам распределения, не правда ли? Проблема в том, что они свалились на нас чудесным образом. Возникает вопрос, можно ли их было не угадать, а вывести как-нибудь? И как быть, если даны не эти конкретные, а какие-то другие ограничения?

Оказывается, что при некоторых не очень обременительных ограничениях ответ можно записать с помощью распределений экспоненциального класса. Давайте же познакомимся с ними поближе.

Экспоненциальное семейство распределений

Говорят, что семейство распределений относится к экспоненциальному классу, если оно может быть представлено в следующем виде:

$$\color{#348FEA}{p(x\vert\theta) = \frac1{h(\theta)}g(x)\cdot\exp\left(\theta^Tu(x)\right)}$$

где $\theta$ — вектор вещественнозначных параметров (различные значения которых дают те или иные распределения из семейства), $h, g > 0$, $u$ — некоторая вектор-функция, и, разумеется, сумма или интеграл по $x$ равняется единице. Последнее, в частности, означает, что

$$h(\theta) = \int g(x)\exp\left(\theta^Tu(x)\right)dx$$

(или сумма в дискретном случае).

Пример 1. Покажем, что нормальное распределение принадлежит экспоненциальному классу. Для этого мы должны представить привычную нам функцию плотности

$$p(x \vert \mu, \sigma^2) = \frac{1}{\sqrt{2\pi}\sigma}\exp\left(-\frac{(x-\mu)^2}{2\sigma^2}\right)$$

в виде

$$p(x\vert\theta) = \frac{g(x)\cdot\exp\left(\sum_i\text{(параметр)}_i\cdot\text{(функция от x)}_i\right)}{\text{что-то, не зависящее от $x$}}$$

Распишем

$$\frac{1}{\sqrt{2\pi}\sigma}\exp\left(-\frac{(x-\mu)^2}{2\sigma^2}\right) = \frac{1}{\sqrt{2\pi}\sigma}\exp\left(-\frac1{2\sigma^2}x^2 + \frac{\mu}{\sigma^2}x + \frac{\mu^2}{2\sigma^2}\right)=$$

$$\frac{\exp\left(-\frac1{2\sigma^2}x^2 + \frac{2\mu}{2\sigma^2}x\right)}{\sqrt{2\pi}\sigma\exp\left(-\frac{\mu^2}{2\sigma^2}\right)}=$$

Определим

$$u_1(x) = x,\qquad u_2(x) = x^2$$

$$\theta_1 = \frac{\mu}{\sigma^2},\quad \theta_2 = -\frac1{2\sigma^2}$$

$$h(\theta) = \sqrt{2\pi}\sigma\exp\left(-\frac{\mu^2}{2\sigma^2}\right)$$

Если теперь всё-таки честно выразить $h$ через $\theta$ (это мы оставляем в качестве лёгкого упражнения), то получится

$$p(x \vert \mu, \sigma^2) = \frac1{h(\theta)}\exp\left(\theta^Tu(x)\right)$$

В данном случае функция $g(x)$ просто равна единице.

Пример 2. Покажем, что распределение Бернулли принадлежит экспоненциальному классу. Для этого попробуем преобразовать функцию вероятности (ниже $x$ принимает значения $0$ или $1$):

$$P(x \vert p) = p^x(1 - p)^{1 - x} = \exp\left(x\log{p} + (1 - x)\log(1 - p)\right)$$

Теперь мы можем положить $u(x) = \left(x, 1 - x\right)$, $\theta = \left(p, 1 - p\right)$, и всё получится. Единственное, что смущает, — это то, что компоненты вектора $u(x)$ линейно зависимы. Хотя это не является формальной проблемой, но всё же хочется с этим что-то сделать. Исправить это можно, если переписать

$$p^x(1 - p)^{1 -x} = (1 - p)\exp\left(x\log{p} + (-x)\log(1 - p)\right) =$$

$$=(1 - p)\exp\left(x\log{\frac{p}{1 - p}}\right)$$

и определить уже минимальное представление с $u(x) = x$, $\theta = \log{\frac{p}{1 - p}}$ — мы ведь уже сталкивались с этим выражением, когда изучали логистическу регрессию, не так ли?

Вопрос на подумать. Принадлежит ли к экспоненциальному классу семейство равномерных распределений на отрезках $U[a, b]$? Казалось бы, да: так как:

$$p(x) = \frac{1}{b - a}\mathbb{I}_{[a,b]}(x)\exp(0)$$

В чём может быть подвох?

Ответ (не открывайте сразу; сначала подумайте сами!)

Нет, не принадлежит. Давайте вспомним, как звучало определение экспоненциального семейства. Возможно, вас удивило, что там было написано не «распределение относится», а «семейство распределений относится».

Это важно: ведь семейство определяется именно различными значениями $\theta$, и если нас интересует семейство равномерных распределений на отрезках, определяемое параметрами $a$ и $b$, то они не могут быть в функции $g(x)$, они должны быть под экспонентой, а экспонента ни от чего не может быть равна индикатору.

При этом странное и не очень полезное семейство с нулём параметров, состоящее из одинокого распределения $U[0,1]$ можно считать относящимся к экспоненциальному классу: ведь для него формула

$$p(x) = \mathbb{I}_{[0,1]}(x)\exp(0)$$

будет работать.

Как мы увидели, к экспоненциальным семействам относятся как непрерывные, так и дискретные распределения. Вообще, к ним относится большая часть распределений, которыми Вам на практике может захотеться описать $Y \vert X$.

В том числе:

• нормальное;
• распределение Пуассона;
• экспоненциальное;
• биномиальное, мультиномиальное (с фиксированным числом испытаний);
• геометрическое;
• $\chi^2$-распределение;
• бета-распределение;
• гамма-распределение;
• распределение Дирихле.

К экспоненциальным семействам не относятся, к примеру:

• равномерное распределение на отрезке;
• $t$-распределение Стьюдента;
• распределение Коши;
• смесь нормальных распределений.

MLE для семейства из экспоненциального класса

Возможно, вас удивил странный и на первый взгляд не очень естественный вид $p(x\vert\theta)$. Но всё не просто так: оказывается, что оценка максимального правдоподобия параметров распределений из экспоненциального класса устроена очень интригующе.

Запишем функцию правдоподобия выборки $X = (x_1,\ldots,x_N)$:

$$p(X\vert\theta) = h(\theta)^{-N}\cdot\left(\prod_{i=1}^Ng(x_i)\right)\cdot\exp\left(\theta^T\left[\sum_{i=1}^Nu(x_i)\right]\right)$$

Её логарифм равен

$$l(X\vert\theta) = -N\log{h(\theta)} + \sum_{i=1}^N\log{g(x_i)} + \theta^T\left[\sum_{i=1}^Nu(x_i)\right]$$

Дифференцируя по $\theta$, получаем

$$\nabla_{\theta}l(X\vert\theta) = -N\nabla_{\theta}\log{h(\theta)} + \left[\sum_{i=1}^Nu(x_i)\right]$$

Тут нам потребуется следующая

Лемма. $\nabla_{\theta}\log{h(\theta)} = \mathbb{E}u(x)$

Доказательство:

Как мы уже отмечали в прошлом пункте:

$$h(\theta) = \int g(x)\exp\left(\theta^Tu(x)\right)dx$$

Следовательно,

$$\nabla_{\theta}\log{h(\theta)} = \frac{\nabla_{\theta}\int g(x)\exp\left(\theta^Tu(x)\right)dx}{\int g(x)\exp\left(\theta^Tu(x)\right)dx} =$$

$$= \frac{\int u(x)g(x)\exp\left(\theta^Tu(x)\right)dx}{h(\theta)} =$$

$$=\int u(x)\cdot\frac1{h(\theta)}g(x)\exp\left(\theta^Tu(x)\right)dx = \mathbb{E}u(x)$$

Кстати, можно ещё доказать, что

$$\frac{\partial}{\partial \theta_i\partial\theta_j}\log{h(\theta)} = \text{Cov}(u_i(x), u_j(x))$$

Приравнивая $\nabla_{\theta}l(X\vert\theta)$ к нулю и применяя лемму, мы получаем, что

$$\color{#348FEA}{\mathbb{E}u(x) = \frac1N\left[\sum_{i=1}^Nu(x_i)\right]}$$

Таким образом, теоретические матожидания всех компонент $u_i(x)$ должны совпадать с их эмпирическими оценками, а метод максимального правдоподобия совпадает с методом моментов для $\mathbb{E}u_i(x)$ в качестве моментов.

И в следующем пункте выяснится, что распределения из семейств, относящихся к экспоненциальному классу, это те самые распределения, которые имеют максимальную энтропию из тех, что имеют заданные моменты $\mathbb{E}u_i(x)$.

**Пример.**Рассмотрим вновь логнормальное распределение:

$$p(x) = \frac1{x\sqrt{2\pi\sigma^2}}\exp\left(-\frac{(\log{x} - \mu)^2}{2\sigma^2}\right) =$$

$$=\frac1{x\sqrt{2\pi\sigma^2}}\exp\left(-\frac1{2\sigma^2}\log^2{x} + \frac{\mu}{\sigma^2}\log{x} - \frac{\mu^2}{2\sigma^2}\right) =$$

$$=\frac1{x\sqrt{2\pi\sigma^2}\exp\left(\frac{\mu^2}{2\sigma^2}\right)}\exp\left(\underbrace{\frac{\mu}{\sigma^2}}_{=\theta_1}\underbrace{\log{x}}_{=u_1(x)} -\underbrace{\frac1{2\sigma^2}}_{=\theta_2}\underbrace{\log^2{x}}_{=u_2(x)} \right) =$$

$$\frac{1}{\sqrt{-\pi\theta_2^{-1}}\cdot\exp{-\frac{\theta_1^2}{4\theta_2}}}\cdot\frac1x\exp\left(\theta_1u_1(x) + \theta_2u_2(x)\right)$$

Как видим, логнормальное распределение тоже из экспоненциального класса. Вас может это удивить: ведь выше мы обсуждали, что для него метод моментов и метод максимального правдоподобия дают разные оценки.

Но никакого подвоха тут нет: мы просто брали не те моменты. В данном случае $u_1(x) = \log{x}$, $u_2(x) = \log^2{x}$, их матожидания и надо брать; тогда для параметров, получаемых из MLE, должно выполняться

$$\mathbb{E}\log{x} = \frac1N\sum_i\log{x_i},\quad \mathbb{E}\log^2{x} = \frac1N\sum_i\log^2{x_i}$$

Матожидания в левых частых мы должны выразить через параметры — и нам для этого совершенно не обязательно что-то интегрировать! В самом деле:

$$\mathbb{E}\log{x} = \frac{\partial}{\partial\theta_1}\log{h(\theta)} =$$

$$=\frac{\partial}{\partial\theta_1}\left(-\frac12\log{\pi} + \frac12\log{\theta_2} - \frac{\theta_1^2}{4\theta_2^2}\right) = -\frac{\theta_1}{2\theta_2^2}$$

$$\mathbb{E}\log^2{x} = \frac{\partial}{\partial\theta_2}\log{h(\theta)} = \frac1{2\theta_2} + \frac{\theta_1^2}{2\theta_2^3}$$

Теорема Купмана-Питмана-Дармуа

Теперь мы наконец готовы сформулировать одно из самых любопытных свойств семейств экспоненциального класса.

В следующей теореме мы опустим некоторые не очень обременительные условия регулярности. Просто считайте, что для хороших дискретных и абсолютно непрерывных распределений, с которыми вы в основном и будете сталкиваться, это так.

Теорема. Пусть $p(x) = \frac1{h(\theta)}\exp\left(\theta^Tu(x)\right)$ — распределение, причём $\theta$ — вектор длины $n$ и $\mathbb{E}u_i(x) = \alpha_i$ для некоторых фиксированных $\alpha_i$, $i=1,\ldots,n$. Тогда распределение $p(x)$ обладает наибольшей энтропией среди распределений с тем же носителем, для которых $\mathbb{E}u_i(x) = \alpha_i$, $i=1,\ldots,n$. При этом оно — единственное с таким свойством: в том смысле, что любое другое распределение, обладающее этим свойством, совпадает с ним почти всюду.

Идея обоснования через оптимизацию.

Мы приведём рассуждение для дискретного случая; в абсолютно непрерывном рассуждения будут по сути теми же, только там придётся дифференцировать не по переменным, а по функциям, и мы решили не ввергать вас в мир вариационного исчисления.

В дискретном случае у нас есть счётное семейство точек $x_1, x_2,\ldots$, и распределение определяется счётным набором вероятностей $p_i$ принимать значение $x_i$. Мы будем решать задачу

$$\begin{cases} -\sum_j p_j\log{p_j}\longrightarrow\max,\\ \sum_jp_ju_i(x_j) = \alpha_i, i = 1,\ldots,n,\\ \sum_jp_j = 1,\\ p_j\geqslant0 \end{cases}$$

Запишем лагранжиан:

$$\mathcal{L} = \sum_j p_j\log{p_j} + \sum_i\theta_i\left(\alpha_i - \sum_jp_ju_i(x_j)\right)+$$

$$+\theta_0\left(\sum_jp_j - 1\right) - \sum_j\lambda_jp_j$$

Продифференцируем его по $p_j$:

$$\frac{\partial\mathcal{L}}{\partial p_j} = \log{p_j} + 1 - \sum_i\theta_iu_i(x_j) + \theta_0 - \lambda_j$$

Приравнивая это к нулю, получаем

$$p_j = \frac{\exp\left(\langle \theta, u(x_j)\rangle\right)}{\exp\left(\lambda_j - \theta_0 - 1\right)}$$

Числитель уже ровно такой, как и должен быть у распределения из экспоненциального класса. Разберёмся со знаменателем.

Легко видеть, что условие $p_j\geqslant0$ заведомо выполнено (ведь тут сплошные экспоненты), так что его можно было выкинуть из постановки задачи оптимизации или, что то же самое, положить $\lambda_j = 0$. Параметр $\theta_0$ находится из условия $\sum_jp_j = 1$, а точнее, выражается через остальные $\theta_i$, что позволяет записать знаменатель в виде $h(\theta)$.

Идея доказательства «в лоб».

Как и следовало ожидать, оно ничем не отличается от того, как мы доказывали максимальность энтропии у равномерного или нормального распределения. Пусть $q(x)$ — ещё одно распределение, для которого

$$\int u_i(x)q(x)dx = \int u_i(x)p(x)dx$$

для всех $i = 1,\ldots,n$. Тогда

$$0\leqslant KL(q\vert\vert p) = \int q(x)\log\left(\frac{q(x)}{p(x)}\right)dx =$$

$$=\underbrace{\int q(x)\log{q(x)}dx}_{-H(q)} - \int q(x)\log{p(x)}dx=$$

$$=-H(q) - \int q(x)\left(-\log{h(\theta)} + \sum_i\theta_iu_i(x)\right)dx =$$

$$=-H(q) - \log{h(\theta)}\underbrace{\int q(x)dx}_{=1=\int p(x)dx} - \sum_i\theta_i\underbrace{\int q(x)u_i(x)dx}_{=\int p(x)u_i(x)dx} =$$

$$=-H(q) - \int p(x)\left(-\log{h(\theta)} + \sum_i\theta_iu_i(x)\right)dx =$$

$$=-H(q) + \int p(x)\log{p(x)}dx = -H(q) + H(p)$$

Таким образом, $H(p)\geqslant H(q)$, причём по уже не раз использованному нами свойству дивергенции Кульбака-Лейблера из $H(p) = H(q)$ будет следовать то, что $p$ и $q$ совпадают почти всюду.

Рассмотрим несколько примеров:

Пример 1. Среди распределений на множестве $\{1,2,3,\ldots\}$ неотрицательных целых чисел с заданным математическим ожиданием $\mu$ найдём распределение с максимальной энтропией.

В данном случае у нас лишь одна функция $u_1(x) = x$, которая соответствует фиксации матожидания $\mathbb{E}x$. Плотность будет вычисляться только в точках $x=k$, $k=1,2,\ldots$ и будет иметь вид

$$p_k = p(k) = \frac1{h(\theta)}\exp\left(\theta k\right)$$

В этой формуле уже безошибочно угадывается геометрическое распределение с $p = 1 - e^{\theta}$. Параметр $p$ можно подобрать из соображений того, что математическое ожидание равно $\mu$. Матожидание геометрического распределения равно $\frac1p$, так что $p = \frac1{\mu}$. Окончательно,

$$p_k = \frac1{\mu}\left(1 - \frac1{\mu}\right)^{k-1}$$

Пример 2. Среди распределений на всей вещественной прямой с заданным математическим ожиданием $\mu$ найдём распределение с максимальной энтропией.

А сможете ли вы его найти? Решение под катом.

Теория говорит нам, что его плотность должна иметь вид

$$p(x) = \frac1{h(\theta)}\exp\left(\theta x\right)$$

но интеграла экспоненты не существует, то есть применение «в лоб» теоремы провалилось. И неспроста: если даже рассмотреть все нормально распределённые случайные величины со средним $\mu$, их энтропии, равные $\frac12 + \frac12\log(2\pi\sigma^2)$ не ограничены сверху, то есть величины с наибольшей энтропией не существует.